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Abstract 

In this paper we analyze the use of time splitting techniques for solving shallow 
water equation. We discuss some properties that these schemes should satisfy so that 
interactions between the source term and the shock waves are controlled. This paper 
shows that these schemes must be well balanced in the meaning expressed by Greenberg 
and Leroux [5]. More specifically, we analyze in what cases it is enough to verify an 
Approximate C-property and in which cases it is required to verify an Exact C-property 
(see PP, [2]). We also include some numerical tests in order to justify our reasoning. 
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1 Introduction 

In this paper, our interest is to analyze time splitting schemes on conservation laws 
with source terms, also called balance laws. A prototype, in one space dimension and under 
certain regularity hypotheses, is given by the following system of partial differential equations 

W(x, t) t + F(W(x, t)) x = G(x, W(x, t)) , (x, t) e R x R+, 
W(x,0) = W (x) , x E R, 

where W : R x IR + — > IR m is the vector of conserved variables, F : M m — > M. m is the vector 
of fluxes and G : M m+1 — > M. m is the source term. 

Recently, there has been some controversy related to the application of time splitting 
techniques on hyperbolic equations involving solutions with shock waves, as is the case of 
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shallow water equations, which are introduced later. Usually, a time splitting numerical 
schemes to solve ([!]) consists of solving consecutively the homogeneous equation 



W{x,t) t + F(W(x,t)) x = 0, 



(2) 



and the ordinary differential equation 



W(x,t) t = G{x,W{x,t)). 



(3) 



LeVeque notices in j9] that such schemes can easily fail by the presence of shock waves 
in solving (T2J). These shock waves involve large changes in the solution which can not be 
captured in solving Q. 

On the other hand, some authors such as Ma, Sun and Yin (see [12]) use a time in- 
tegrating scheme with two-step predictor-corrector sequence quite successfully Striba use 
splitting techniques (see [13]) in meteorology models on the term that represent de Coriolis 
acceleration. In addition, Wicker and Skamarock (see [16]) use time-splitting methods for 
integrating the elastic equations. 

In 1994, Bermudez and Vazquez (see [T], [2]) introduce the concept of Exact C-property 
and Approximate C-property in order to identify numerical schemes with an acceptable level 
of accuracy in the resolution of shallow water equations (well-balanced scheme). Going 
deeper into well-balanced scheme idea, we can find the work of Greenberg and Leroux [5], in 
which they propose a numerical scheme that preserves a balance between the source terms 
and internal forces due to the presence of shock waves. 

Another point of view, more recently, it is provided by Lubich (see [H]), who gives 
an error analysis of Strang-type splitting integrators for nonlinear Schrodinger equations. 
Holdahl, Holden and Lie use an adaptive grid refinement and a shock tracking technique 
to construct a front-tracking method for hyperbolic conservation laws. They combine the 
operator splitting to study the shallow water equations (see [6]). Holden, Karlsen, Risebro 
and Tao (see [8]) show that the Godunov and Strang splitting methods converge with the 
expected rates if the initial data are sufficiently regular. Finally, in this way the reader can 
find a deep study of splitting methods for partial differential equations in [7], where some 
analysis of conservation and balance laws are included. 

We base our analysis on the ideas presented in p], [2] and [5]. From these studies, it 
follows that the numerical scheme used in solving fl2]) and the numerical scheme used in 
solving ([3]) cannot be whatever, even though they have a high degree of accuracy. These 
must be balanced so that interactions between the source term in ([3]) and the shock waves in 
([2]) are controlled. In this framework, we will analyze conditions to be verified by splitting 
schemes in order to avoid spurious oscillations, which are created in this type of equations. 
More specifically, in which cases it is enough to verify an Approximate C-property and in 
which cases it is required to verify an Exact C-property. 

The rest of this paper is organized as follows: in Section 2, we introduce the governing 
equations for the one dimensional shallow water model. We also analyze two kinds of time 
splitting schemes to identify which conditions must satisfy a scheme for solving this type 
of equations. In Section 3, we show the four test problems that we will use to test the 
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performance of the schemes described above. In Section 4, we present the numerical results 
we have obtained. Finally, in Section 5 we reason the final conclusions we have reached by 
using the results obtained with the presented schemes. 

2 The one dimensional shallow water equations 

2.1 Governing equations 

In this section we consider Eq. (CEJ) with 
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where the unknowns of the problem are: the water height is h and the flow per unit length 
is q = hu. Here u is the average vertical speed in the direction of the axis x (see Fig.l). F is 
the flux of conservative variables and g = 9.81ms -2 is the acceleration due to gravity. The 
source term G models the bottom variation given by the function b(x). 




Fig. 1: Shallow water variables. 



Let us consider numerical solvers based upon the decomposition F(W) X = A(W) W x , 
where 

/ 1 \ 

A(W) = n 2 is the Jacobian matrix of F(W). (5) 



hi 
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So that, system (pQ)-(jl]) is hyperbolic {h > 0), then A = XAX 1 , where 



where Ai = \ + \ fgh and A2 = % — \fgh- 
h h 



Our analysis needs the definition of some conservation properties given by Bermudez and 
Vazquez in j2]. Since only the source term involves the bed slope, an inadequate choice of 
the numerical schemes can give spurious oscillations. So, these conservation properties try to 
identify which kind of schemes have good behavior in equations with source term. Bermudez 
and Vazquez characterize the good behavior of the numerical scheme in the manner in which 
the scheme approximates a steady solution representing the state of water at rest. They 
introduce the stationary problem {Problem SP) given by q(x, t) = and h(x, t) = H(x) and 
define the following conservation properties: 

Definition 1 Exact C-property. We say that a scheme satisfies the Exact C-Property 
if it is exact when applied to the Problem SP. 

Definition 2 Approximate C-property. We say that a scheme satisfies the Approxi- 
mate C-Property if it is accurate to the order B(Ax 2 ) when applied to Problem SP. 

When a numerical scheme does not satisfy any of these conservation properties then the 
propagation of spurious oscillations is also present in non stationary problems. 



2.2 Central numerical schemes 



The differential formulation of the homogeneous equation given in ([2]) does not admit dis- 
continuous solutions. So, since these solution are physically relevant in this context, we need 
a suitable formulation of the problem to support discontinuous solutions. In this sense, it is 
usual to consider the following integral formulation 



/ 



(Wdx - F(W)dt) = 0. (7) 

The numerical schemes use usually (0) in order to approximate (j2J). To do that, it is 
introduced a control volume in the space (x, t) of dimensions Ax x At. Next, it is evaluated 
the integral (JTJ) in this volume control 



r i-l/2 



(W(x, t n+1 ) - W(x, t n ))dx + / (F(W(x j+1/2 ,t)) - F{W{x^ 1/2 , t)))dt = 0. 
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Dividing by Ax we obtain 



i r x j+i/2 i r x j+i/2 

W(x,t n )dx 



Ax 



/ W(xX +1 )dx = — / 



Ax 



Thus, we deduce the conservation formula 



W™ +1 — — ^~~[-*j+i/2 — Fj-iM, (8) 



where W a is an average 

1 r x j+i/2 



Ax 

at time t = t n inside the interval 



W] =— / W(x,t n )dx 



whose length is 

Ax = Xj+i/2 - Xj-i/2- 

The flux in (jSJ) can be interpreted as the average in time of the physical flux, i.e., 

1 r +i 

Fj+1/2 ~AtJ n F (W(x j+1/2 ,t))dt. (9) 

Conservative numerical methods for (j2J) are based in (jSJ, and they are determined by the 
expression of the numerical flux Fj +1 / 2 . 

The time splitting numerical schemes (see [H]), in each time step, act as follows: taking 
into account the initial condition W n , we solve d2J) and we obtain W n+1 . Then, by using the 
initial condition W n+1 , we solve fl3]) and obtain W n+l . 

More specifically, if we use the operators 

• A(W n ) = W n+1 such as 

f W t + F(W) x = 0, (x,t)eRx[t n ,t n+1 ], 
[ W(x,t n ) = W n , x G E. 

• = W n+1 such as 

W t = G(x,W), (x,t) GRx [r,r +1 ] 
jy n ->■ < _ !> iy n+1 

W(x,t n ) = W n+1 , x G R. 



Thus, 



In the first step of the splitting procedure, we solve (J2J) at each time step. To do that 
and for the sake of simplicity, we consider the Q-scheme of van Leer (see [J5] ) which uses a 
matrix Q satisfying some properties and the numerical fluxes 



(10) 



where the numerical flux function <fi is given by 



m v) = F{U) | Fiy) - \\Q(u, v)\(v - u). 



A possible choice of the matrix Q can be the Jacobian (jSJ) of the system 
the arithmetic mean, i.e. 

Q(U,V) = A l 



So that, we have 

Qiw; ±1 



X(W? ±1 , Wf) |A(W7 ±1 , W?) I x-\w? ±x , W?) 



evaluated at 
(12) 

(13) 



where X(W? ±1 , W?), \A(W? ±l7 W?)\ and X^iW^ Wf) are evaluated at - 



So, we obtain the numerical fluxes 



^•±1/2 



F(W? ±1 ) + F(W} 



q(w; ±l ) ( T wf±vr; ±1 ). 



(14) 



The second step of the procedure is to solve (j3J). To do that, we get the solutions of 
the homogeneous equation (j2J): h(xj,t n+1 ) and g(xj,t n+1 ) and we solve the following initial 
value ODE problem for each Xj 



d ( h(xj,t) 
dt \ q ( Xj , t ) 

h{xj,t n ) 
q(xj,t n ) 





—gh(xj, t)b'(xj 

h{xj,t n+1 ) 
q(xj,t n+1 ) 



t G [t n ,t 



(15) 



The first equation has the solution: 



h( Xj ,t n+1 ) = h(xj,t n+1 ] 



(16) 
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To solve the second one, we calculate 



dq(xj, t) 
dt 



-gh(xj,t)b'(xj), 



q{xj,t) = - gh(xj,s)b'(xj) ds; t G (t n ,t n+1 ]. 



Therefore, 



qfe, t n+1 ) = q( Xj , t n+1 ) - g ■ b'{xj) ■ / h{ Xj , s) ds. 



Finally, by using a trapezoidal rule, we obtain 



q{x 3) t n+1 ) = q(x v t n+1 ) - g ■ b'( Xj ) ■ ^(h(x v t n ) + h(x 3 ,t n+1 )). (17) 



Remark 1 We denote the time splitting scheme given by (TJ\)- [Tfy) -(T7\) as Q-tral. 
Proposition 1 The numerical scheme Q-tral satisfies an approximate C-property. 
Proof 



Problem SP satisfies Ai = ygh and A2 = —ygh. In addition (see Fig. 1), b(x) + h(x,t) ■ 
b(x) + H(x) = A. Therefore, 

b( Xj ) = A - H( Xj ). (li 



In order to solve (jzj) with ffl4l) . we compute 



Q(W, 



3+1 



Then, 



|Ax|A 2 - |A 2 |A! -lAil + IAaj 



( . I k 3+l + 5 



A 2 - A x V AiA 2 (|Ai| - |A 2 |) -AxlAxI + A 2 |A 2 



hjjf-\ + hj 



~2 / 



/ 



(h j+1 - hj) \ g 



h j+1 + hj \ 



Analogously, we obtain 



■> 2 



(hj - h^) \J~g 



hj + hj-i \ 



5 + h U 



Finally, using (J5J), we have the numerical approximation of Eq. fl2]) 



W, 



n+l 





1- 







At 
2 Ax 



hj + hj-i 



V 



gAt_ 
4 Ax 



W + i ~ hU) 



) 



Next, to have the numerical approximation of Eq. fl3]), we use f fl6|) and (fT7j) . By substi- 
tuting, we obtain 



At 
2 Ax 



hj + hj-i 



(19) 



and 



From (fl9l) . we deduce that = foj is accurate to the order B(Ax 2 ). And from (T50 



and (I18p we have that g"' = is exact when we take 



2 Ax 



Thus, the proof ends. 



□ 



Remark 2 If we use another convergent quadrature rule to approximate the integral / h(xj, 

Jt n 

then Proposition 1 is also satisfied. 



2.3 Upwind numerical schemes 

It is well-known that conservation laws with source terms can be solved with high accuracy 
and discontinuities are well captured by using upwind schemes. For instance, you can see the 
works of LeVeque and Yee [JO] or Vazquez- Cendon [15]. Furthermore, this kind of schemes 
have bigger stability regions than schemes using centered approximations of the source term 
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(see [2]). In this section, we will use some approaches in order to build an upwind scheme 
using some ideas contained in [15] and making changes to solve fl3]). 

In order to solve ([3]), we need to introduce two numerical source functions, Gl on the 
left and G R on the right, to upwind the source term. We also use two matrices D L and D R 
with the aim of making clear the upwind process. Below we explain the process in detail. 

We will use, in each time step, the average of solutions of the following equations: 



w t = G L (x,w), (x,t) six [r,r +1 ] 

W{x,t n ) = W n+ \ xeR. 



(21) 



and 



w t = g r (x, w), {x,t) gMx [r,r +1 ] 

W(x,t n ) = W n+ \ 168, 



(22) 



where 



G L (x, W) = D L {W L ) G(x, W) and G R (x, W) = D R {W R ) G(x, W)\ 



and functions W l and W R change in each computacional cell in the discretization problem 

— Wj + W 3 . x — W j+1 + Wj 

W L = — — and W R = -. 

2 2 

Thus, we obtain the solution of fT2T]) -f l22|) . W£ +1 and W R +1 respectively. Then, in each 
time step, we take as solution of (J3J) the following average 



W 



n+l 



W n+l + W n+l 



In other words, for the discretization of the problem, we can substitute the ODE in ffT5"j) 



by 



d_ 

dt 



h(xj, t) 
q(xj,t) 



I),.(\Y l .;)G(x l .\V J 



^21 *^22 



-gh(xj,t)b'(xj 



-d^ghM 



Therefore, instead of ( fl5l) . we have two ODE 



d_ 

dt 



h(xj,t r 
q(xj,t n 



-d^ghjb'j 
-d^ghjb'j 

" /i(x j ,r +1 ) 

q(xj,t n+1 ) 



t e [t n ,t 



n j.n+1] 



(23) 
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and analogously from (|22|) we obtain 



d_ 

dt 



h{xj,t n ) 
q(xj,t n ) 



-d? 2 gh j b' j 
-dgghjb'j 

h(xj,t n+1 ) 
fe,t n+1 ) 



t e [t n ,t 



n j.n+11 



(24) 



where d{\ and d 22 are coefficients corresponding to matrix Dr. 
Finally, we compute the solutions of (125]) and (123]) , W^ 1 and W^ 1 . 

Remark 3 Matrices Dl and Dr must be chosen in coordination with the numerical method 
used in solving |1P to get a well balanced scheme. For example, regarding ( fl^l ), we can take 



D L =(I+ \Q\Q- 1 ) and D R =(I- \Q\Q~ 1 ) 
In this way we obtain consistency, since Dl + Dr = 21, we have 

G L + G R D l G + D R G (D l + D R )G 



(25) 



G. 



2 2 2 

Remark 4 We denote the time splitting scheme given by (T^- f^) -(24\)- f2l)\) as Q-tra2. 

Proposition 2 The numerical scheme Q-tra2 satisfies an exact C-property. 

Proof 

Acting in a similar way to Proposition 1, we have 



W. 



( 

I hi \ 

n+l 

v 

On the other hand, 
1 



At 
^ + 2A^ 



D x {W l 



Now, taking in ( 123 



(h j+1 - hj)\l g hj+1 ^ hj - (hj - h^) 



hj + hj-i\ ^ 



gAt_ 
4 Ax 



W + i ~ hU) 



and D 2 (W Rj 



h j+1 +hj 

y 2 



b'{Xj) = h M Z^zO and h(Xj,t) = h > 
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to integrate, we have 



and 



h Ll = h 3 - —(bi - bj-JJg 3 2 3 (26) 



QLj = Qj ~ ^^9(bj - bj^ihj + hj-x) (27) 



Acting in a similar way to ( 124)) . we have 



h R] =h J -^(b J+1 -b 3 ) ] j 9 f ^±^ (28) 



and 

At 



QRj = Qj ~ ~ b j)( h J+i + h i) ( 29 ) 



In addition, from (fT8|) we have bj = A — hj, then 



+1 = h Lj + h Rj = ^ = g Lj + g Rj = q 

7 2 J 2 

This concludes the proof. 

□ 



3 Numerical experiments 



In this section, we will discuss four test problems for shallow water equations. The first test 
problem is the dam-break problem with variable topography, which is one of the most basic 
problems with source term. The second one is a stationary problem, also with source term 
that represents a smooth bottom. The third one represents a tidal wave flow propagating 
over an irregular topography. This test was discussed by Bermudez and Vazquez [2] and it 
is one of the most popular test to check the performance of a numerical scheme which tries 
to be effective in solving shallow water equations. It represents a severe test regarding to 
irregular topography and the long time over which is applied. Finally, the last test represents 
a numerical simulation of a tidal wave on the shoreline with friction effects, which introduces 
a wet /dry front. 
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Test 1: Dam-break with variable topography 



We consider the one dimensional shallow water equations fll]) on the domain (x,t) £ [0, 1] x 
R + . The initial data are 



1 / / 1 
- cos IOtt x 

b(x) = \ 8 V V 2 
0, otherwise; 



1 2 3 
1, - < x < -, 

5 5' W(x,0) 



1 - b{x) 


0.5 - b(x) 




, x < „ 



2 ■ 



Fig. 2 shows numerical results of Test. 1 for schemes Q-tral and Q-tra2. In this test, we 
have chosen cfl = 0.5 and t = 0.5 on 200 computational cells. The number of cells is large 
and the value of t is small, so under these conditions, one can see that the performance of 
both schemes is similar. 



Initial condition 



"W 

free surface 


rail 

free surface 


/ 


\b(x) 



t = 0.5 



oooooo Q-tr8/2 
- Q-tral 



0.2 0.4 0.6 0.8 



0.2 0.4 0.6 0.8 



Fig. 2: Test 1. Initial condition and numerical result for t = 0.5. 



Test 2: Stationary problem with smooth bottom 

We consider also Eq. (j4j) on the domain (x, t) £ [0, 1] x M + . The bottom function b(x) is the 
same of test 1 and the initial data are 



1 - b(x) 

W(x,0) = ( q | , <x < 1. 



This problem has a stationary solution W(x, t) = W(x, 0), < x < 1 and t > 0. 
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In this test, we have chosen cfl = 0.5 and t = 0.25 on 50 computational cells, it shows 
spurious oscillations for Q-tral (see Fig. 3). However, Q-tra2 scheme has a good behavior. 



o°o 
o o 
o o 



- Analytical surface 



oooooo 
OOOOOO 



Q-tral 
Q-tra2 



0.2 0.4 0.6 0.6 



Fig. 3: Test 2. Numerical result for t = 0.25. 



Test 3: Tidal wave flow with irregular topography 

We consider Eq. fll]) on the domain (x,t) G [0, 1500] x [0, 10800], where length is in meters 
and time in seconds. The bottom function b(x) is shown in Fig. 4, and the initial data are 



H(x) 

W(x, 0) = [ ) , H(x) = H(0) - b(x), 0<x< 1500 and H(0) = 16. 



Moreover, we have the following boundary conditions 



h(Q,t) = H{0) + <p{t), 



ip(t) = 4 + 4 sin ( 7r 
g(1500,t) = 0, 



At 



86400 



} < t < 10800. 



Function <p(i) simulates a tidal wave of 4m amplitud. 
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For this problem, we can obtain an asymptotic analytical solution (see [2]) given by 
h(x,t) = H(x) + <p(t) and it satisfies h(x, 10800) = 20, < x < 1500. 

In order to check the behaviour of Q-tra2 scheme, we have chosen cfl = 0.9 on 100 
computational cells. Test. 3 is an strong test, it uses a long time t = 10800 seconds and a 
complex topography, however an excellent performance of this scheme is demonstrated (see 
Fig. 4). 




Fig. 4: Test 3. Tidal wave flow over an irregular topography using Q-tra2. 



Test 4: Tidal wave on the shoreline with friction effects 



In this test we take in account the bottom friction. Manning law is used to model friction 
between fluid and bottom, so a new element appears in the source term. As a consequence, 
the source term is given by 







G(x,W) 



h-l 



-ghb'(x) - gqM 2 

where M is the Manning coefficient. 

We consider the domain (x, t) e [0, 6] x R + and the bottom function 

0.00125x + 0.0125, < x < 3, 
0.162(3 - 3) + 0.01625, 3 < x < 6. 



b(x) 



Here, length is in meters, time in seconds and the reference level is located at 40 cm. 
The tidal wave is simulated by introducing a discharge q = 0.8m 3 /s in the boundary x = 
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during the first 0.2 seconds, from that time a vertical wall condition is introduced. In the 
boundary x = 6 is also introduced a vertical wall condition to preserve the mass conservation 
in the domain. 



On the other hand, this test involves a wet/dry front, which complicates the stability of 
the scheme. Specific treatment is required to remove the spurious oscillations, the reader 
can find all detailed in j3] and pE]. For brevity, we summarize here two main actions of this 
treatment: 

1.- Redefinition of the discretized bottom function to avoid the appearance of spurious 
pressure forces. 

• If Ij is dry, Ij_i wet and hj_\ + bj_x < hj + bj, then bj = + foj-i- 

• If Ij is wet, dry and hj_\ + bj_\ > hj + bj, then bj = — hj. 

2 - Simulate the fact that the discharge is zero when the fluid across a wet/dry front. 

• If Ij is dry, then q] +1 = 0. 

• If Ij is wet, estimate of < and Jj_i is dry, then qj +1 = 0. 

• If Ij is wet, estimate of > and Ij + \ is dry, then qj +1 = 0. 



In addition, friction causes spurious oscillations in g". However splitting schemes can 
avoid these effects by using a semi-implicit discretization. Below, we summarize the details. 



In order to compute W^ 1 , we consider Eq. (12"B|) with friction 



d_ ( h j 

dt { g. 



h(xj,t n ) 

\ q(xj,t n ) 



-d\ 2 ghjb'j - d\ 2 gqjM 2 



-d\ 2 ghjb'j - d{ 2 gqjM 2 



h(xj,t n+1 ) 



( 1l 

h\ 



{hi, 
(hi, 



4 \ 



; te[t n ,t 



(30) 



q(xj,t 



J^+l^ 



Solving, we have 



h Lj = hj - d\ 2 Atg 



ks+ n hs - 1 V i -dlAtg( hi + hi - 1 



* + V- 1 M 2 



Qj + gj-i 
hj + hj-i 



(31) 
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and 



qL 3 = qj - dl 2 Atg kj + J 1 b' j 



dl 2 Atg 



hj + hj. 



Qj + Qj-i 



(32) 



Now, from Eq. (1321) we obtain that 



g, - d\ 2 Atg^±^b> 3 



1 + d\ 2 Atg 



M 



qj + ^ 



hj + hj_i 



(33) 



On the other hand, in order to compute W^ 1 , we act analogously and obtain 



h 



hj - dl 2 Atg^±^b' - d\ 2 Atg ( hj+l + ^ 



and 



QRi 



gj+i + Qj M 2 



Qj ~ d\ 2 Atg- 



ij+i 



Qj+i + Qj 



hj+i + hj 



d\ 2 Atg [ h±l±k V* A 



qj+i + qj 



hj+i + hj 



(34) 



(35) 



Remark 5 VKe denote the time splitting scheme given by p4\)- f3l\ )- l[^) -(34^ and < f^3)j as 
Q-tra3. 



Test. 4 results using Q-tra3 scheme are shown in Fig. 5. This test has all the ingredients 
of a real problem in shallow water: irregular topography, friction effects between bottom 
and water, and advancing front wet/dry. In Fig. 5, we can see the progress of the tidal wave 
for different time values: for t — 1 the wave goes to head, in t = 2 the wave reaches the 
shoreline, in t = 3 the wave hits the wall and in t — 4, it returns toward x — 0. In t — 5, we 
can see how the wave returns back with negative velocity and it crashes with the wave from 
x = with positive velocity. However, in front of this critical situation, the behavior of the 
scheme is stable. 

In this test, we have used a Manning coefficient M = 0.015 and cfl = 0.5 on 250 
computational cells. It is appropriate to note that if you do not use a semi-implicit scheme, 
we obtain spurious oscillations and overflow computations. 
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t = 1 



t = 2 




4 Concluding remarks 

Some upwind time splitting schemes have been developed for solving 1-D shallow water 
equations. In some test problems, our analysis shows that it is not enough to verify an 
approximate C-property to obtain good approximations; we need that an exact C-property 
to be verified. Perhaps, for some undemanding tests, it is not necessary, but when the 
test is demanding, then schemes that only satisfy an approximate C-property can bring to 
numerical instabilities. 

We have shown that time splitting schemes for solving shallow water equations must be 
well balanced in the following sense: the scheme used to solve the homogeneous equation 
Eq. 02]) and the scheme used to integrate the time-dependent ODE Eq. ((3]) cannot be 
whatever. They must be linked so that an exact C-property is satisfied. Apparently, Q-tra2 
is computationally very expensive. But this is not true, since carried out calculations to 
obtain the numerical approximation of Eq. ([2]) are used to obtain Dl and Dr matrices. 
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Moreover, time splitting schemes could give us the possibility to choose a semi-implicit 
or implicit solver in order to obtain stability when Eq. fl3]) is integrated. In Test. 4, we have 
used a semi-implicit scheme that allowed us to obtain stability when shallow water model 
takes into account friction effects between fluid and bottom. 

In summary, we have presented three time splitting schemes with different properties, 
the second of which is the best suited to treat severe test on shallow water equations. The 
semi-implicit Q-tra3 scheme shows the versatility that can provide splitting techniques for 
solving challenges of stiffness that present some partial differential equations. 
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